Morphological divergence in fragmented populations: No evidence of fitness benefit

Study Design

Setup

Load all packages

# libraries
library(readr)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(gganimate)
library(gifski)
library(dplyr)
library(tidyr)
library(ggthemes)
library(ggpubr)
library(lme4)
library(ggeffects)
library(lmerTest)
library(car)
library(ggplot2)
library(survival)
library(coxme)
library(survminer)
library(ggfortify)
library(broom)
library(ehahelper)
library(dena)
library(tidyr)

sessionInfo()

sessionInfo()

R version 4.3.0 (2023-04-21)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/London
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] dena_0.1.0         ehahelper_0.3.9999 broom_1.0.4        ggfortify_0.4.16  
 [5] survminer_0.4.9    coxme_2.2-18.1     survival_3.5-5     car_3.1-2         
 [9] carData_3.0-5      ggeffects_1.2.2    ggpubr_0.6.0       ggthemes_4.2.4    
[13] tidyr_1.3.0        dplyr_1.1.2        gifski_1.6.6-1     gganimate_1.0.8   
[17] ggplot2_3.4.2      sjlabelled_1.2.0   sjmisc_2.8.9       sjPlot_2.8.14     
[21] readr_2.1.4        bdsmatrix_1.3-6    lmerTest_3.1-3     lme4_1.1-33       
[25] Matrix_1.5-4      

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.0    farver_2.1.1        fastmap_1.1.1       tweenr_2.0.2       
 [5] bayestestR_0.13.1   sjstats_0.18.2      digest_0.6.31       estimability_1.4.1 
 [9] lifecycle_1.0.3     magrittr_2.0.3      compiler_4.3.0      rlang_1.1.1        
[13] sass_0.4.6          progress_1.2.2      tools_4.3.0         utf8_1.2.3         
[17] yaml_2.3.7          eha_2.10.3          data.table_1.14.8   knitr_1.42         
[21] ggsignif_0.6.4      prettyunits_1.1.1   abind_1.4-5         withr_2.5.0        
[25] purrr_1.0.1         numDeriv_2016.8-1.1 grid_4.3.0          fansi_1.0.4        
[29] xtable_1.8-4        colorspace_2.1-0    emmeans_1.8.5       scales_1.2.1       
[33] MASS_7.3-60         insight_0.19.1      cli_3.6.1           mvtnorm_1.1-3      
[37] rmarkdown_2.21      crayon_1.5.2        generics_0.1.3      km.ci_0.5-6        
[41] rstudioapi_0.14     performance_0.10.3  modelr_0.1.11       tzdb_0.3.0         
[45] minqa_1.2.5         cachem_1.0.8        stringr_1.5.0       splines_4.3.0      
[49] survMisc_0.5.6      vctrs_0.6.2         boot_1.3-28.1       jsonlite_1.8.4     
[53] hms_1.1.3           rstatix_0.7.2       jquerylib_0.1.4     glue_1.6.2         
[57] nloptr_2.0.3        cowplot_1.1.1       stringi_1.7.12      gtable_0.3.3       
[61] munsell_0.5.0       tibble_3.2.1        pillar_1.9.0        htmltools_0.5.5    
[65] KMsurv_0.1-5        R6_2.5.1            evaluate_0.20       lattice_0.21-8     
[69] backports_1.4.1     bslib_0.4.2         Rcpp_1.0.10         gridExtra_2.3      
[73] nlme_3.1-162        xfun_0.39           zoo_1.8-12          pkgconfig_2.0.3    

Data files

Data file ADD_I is available on the github page.

Morphometric and Island Comparison

Preparing Data

# read the data file "ADD_I"
ADD_I <- read_csv("Tables/ADD_I.csv", col_types = cols(Island = col_factor(levels = c("CN","AR", "CE", "DS", "FR"))))
# create the earliest born individuals on each island
Isfy <- ADD_I %>% select(founderyear,Island) %>% filter(!duplicated(Island))
# create a data frame with means and sd of each morphometric
ADD_mean <- ADD_I %>% group_by(Island,BirthYear) %>% summarise(RTsd = sd(RightTarsus,na.rm=TRUE),WLsd = sd(WingLength,na.rm=TRUE),BMsd = sd(BodyMass,na.rm=TRUE),HBsd=sd(HeadBill,na.rm=TRUE),RightTarsus=mean(RightTarsus,na.rm=TRUE),WingLength=mean(WingLength,na.rm=TRUE),BodyMass=mean(BodyMass,na.rm=TRUE),HeadBill=mean(HeadBill,na.rm=TRUE))

Linear Mixed Effect Models

Tarsus Length

####RightTarsus ----
IRTlmer <- lmer(RightTarsus ~  Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IRTlmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: RightTarsus ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +  
    (1 | BirdID)
   Data: ADD_I

REML criterion at convergence: 10511

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.0560 -0.4445  0.0052  0.4324  5.5726 

Random effects:
 Groups   Name        Variance  Std.Dev.
 BirdID   (Intercept) 0.3048654 0.55215 
 Observer (Intercept) 0.0562847 0.23724 
 Ageclass (Intercept) 0.0009875 0.03142 
 Residual             0.1407944 0.37523 
Number of obs: 6154, groups:  BirdID, 3258; Observer, 82; Ageclass, 5

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     2.417e+01  6.083e-02  2.713e+02 397.385  < 2e-16 ***
IslandAR       -7.215e-02  1.686e-01  3.839e+03  -0.428 0.668783    
IslandCE       -2.510e-01  9.829e-02  3.925e+03  -2.553 0.010704 *  
IslandDS        2.519e-01  8.080e-02  3.408e+03   3.117 0.001842 ** 
IslandFR        1.869e-01  1.236e-01  2.500e+03   1.512 0.130545    
fYear           8.838e-03  2.530e-03  2.174e+03   3.493 0.000488 ***
SexMale         1.596e+00  2.244e-02  3.211e+03  71.130  < 2e-16 ***
IslandAR:fYear -6.760e-04  1.024e-02  3.539e+03  -0.066 0.947388    
IslandCE:fYear  1.764e-03  6.196e-03  4.079e+03   0.285 0.775832    
IslandDS:fYear -2.910e-03  5.840e-03  4.085e+03  -0.498 0.618312    
IslandFR:fYear  1.554e-03  1.756e-02  1.981e+03   0.089 0.929486    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) IslnAR IslnCE IslnDS IslnFR fYear  SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR    -0.163                                                               
IslandCE    -0.309  0.111                                                        
IslandDS    -0.473  0.113  0.218                                                 
IslandFR    -0.331  0.081  0.150  0.225                                          
fYear       -0.779  0.188  0.339  0.504  0.366                                   
SexMale     -0.177 -0.029  0.008 -0.033 -0.024 -0.016                            
IslndAR:fYr  0.149 -0.922 -0.096 -0.101 -0.076 -0.200  0.026                     
IslndCE:fYr  0.261 -0.094 -0.915 -0.173 -0.126 -0.323 -0.010  0.097              
IslndDS:fYr  0.314 -0.065 -0.138 -0.849 -0.208 -0.379  0.008  0.069  0.123       
IslndFR:fYr  0.102 -0.025 -0.050 -0.108 -0.850 -0.143  0.020  0.029  0.048  0.215

#predicting the model 
IRTlmerdata <- ggpredict(IRTlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,RightTarsus=predicted,Island=group,Sex=facet)
IRTlmerdata2 <- merge(IRTlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear) %>% summarise(RightTarsus = mean(RightTarsus, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))

#plotting the model (shown as ggarrange below)
IRTmod <- ggplot(IRTlmerdata2, aes(x = BirthYear, y = RightTarsus, fill = Island)) + 
  geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=RightTarsus,colour=Island), size = 2) +
  geom_errorbar(data=ADD_mean,aes(ymin=RightTarsus-RTsd,ymax=RightTarsus+RTsd, colour=Island),alpha=0.7) +
  geom_smooth(aes(colour = Island),method = "loess", se = FALSE) + 
  geom_ribbon(aes(ymin=RightTarsus-std.error,ymax=RightTarsus+std.error),alpha=0.5) +
  labs(x = "Birth Year", y= "Tarsus Length (mm)") +
  xlim(1990,2022) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind() + scale_fill_colorblind()

Body Mass

#### Body Mass ----
IBMlmer <- lmer(BodyMass ~  Island*fYear + Sex  + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IBMlmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: BodyMass ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +  
    (1 | BirdID)
   Data: ADD_I

REML criterion at convergence: 16677.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6968 -0.5705 -0.0271  0.5440  4.2625 

Random effects:
 Groups   Name        Variance Std.Dev.
 BirdID   (Intercept) 0.24508  0.4951  
 Observer (Intercept) 0.05818  0.2412  
 Ageclass (Intercept) 0.14657  0.3828  
 Residual             0.53384  0.7306  
Number of obs: 6554, groups:  BirdID, 3398; Observer, 91; Ageclass, 5

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     1.438e+01  1.846e-01  5.339e+00  77.909 2.29e-09 ***
IslandAR        1.189e+00  1.570e-01  3.120e+03   7.577 4.63e-14 ***
IslandCE        2.792e-01  1.214e-01  3.666e+03   2.299  0.02158 *  
IslandDS        7.385e-01  1.037e-01  1.096e+03   7.122 1.92e-12 ***
IslandFR        5.264e-01  1.668e-01  1.471e+03   3.156  0.00163 ** 
fYear          -2.623e-03  3.107e-03  5.626e+02  -0.844  0.39893    
SexMale         1.433e+00  2.635e-02  2.907e+03  54.405  < 2e-16 ***
IslandAR:fYear -7.870e-02  1.043e-02  2.346e+03  -7.546 6.38e-14 ***
IslandCE:fYear -3.708e-02  7.935e-03  4.479e+03  -4.673 3.05e-06 ***
IslandDS:fYear -5.069e-02  7.678e-03  2.097e+03  -6.601 5.14e-11 ***
IslandFR:fYear -5.400e-02  2.379e-02  1.097e+03  -2.270  0.02343 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) IslnAR IslnCE IslnDS IslnFR fYear  SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR    -0.071                                                               
IslandCE    -0.107  0.098                                                        
IslandDS    -0.176  0.107  0.184                                                 
IslandFR    -0.113  0.071  0.120  0.187                                          
fYear       -0.307  0.188  0.294  0.462  0.320                                   
SexMale     -0.068 -0.016  0.007 -0.035 -0.027 -0.017                            
IslndAR:fYr  0.065 -0.869 -0.082 -0.095 -0.067 -0.205  0.016                     
IslndCE:fYr  0.092 -0.086 -0.910 -0.144 -0.097 -0.283 -0.007  0.092              
IslndDS:fYr  0.117 -0.064 -0.114 -0.847 -0.181 -0.345  0.015  0.069  0.098       
IslndFR:fYr  0.035 -0.024 -0.039 -0.097 -0.868 -0.127  0.023  0.027  0.035  0.201

#predicting the model
IBMlmerdata <- ggpredict(IBMlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,BodyMass=predicted,Island=group,Sex=facet)
IBMlmerdata2 <- merge(IBMlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear)  %>% group_by(Island,BirthYear) %>% summarise(BodyMass = mean(BodyMass, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))

#plotting the model
IBMmod <- ggplot(IBMlmerdata2, aes(x = BirthYear, y = BodyMass, fill = Island)) + 
  geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=BodyMass,colour=Island), size = 2) +
  geom_errorbar(data=ADD_mean,aes(ymin=BodyMass-BMsd,ymax=BodyMass+BMsd, colour=Island),alpha=0.7) +
  geom_smooth(aes(colour = Island),method = "lm", se = FALSE) + 
  geom_ribbon(aes(ymin=BodyMass-std.error,ymax=BodyMass+std.error),alpha=0.5) +
  labs(x = "Birth Year", y= "Body Mass (g)") +
  xlim(1990,2022) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind() + scale_fill_colorblind()

Wing Length

#### Wing Length ----
IWLlmer <- lmer(WingLength ~  Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IWLlmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: WingLength ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +  
    (1 | BirdID)
   Data: ADD_I

REML criterion at convergence: 23050.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8213 -0.5350  0.0161  0.5335  4.0139 

Random effects:
 Groups   Name        Variance Std.Dev.
 BirdID   (Intercept) 0.95940  0.9795  
 Observer (Intercept) 0.40157  0.6337  
 Ageclass (Intercept) 0.08901  0.2983  
 Residual             1.26489  1.1247  
Number of obs: 6510, groups:  BirdID, 3408; Observer, 91; Ageclass, 5

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     6.592e+01  1.905e-01  1.616e+01 346.012  < 2e-16 ***
IslandAR        2.490e-01  2.795e-01  3.717e+03   0.891  0.37296    
IslandCE       -1.869e-01  2.098e-01  4.116e+03  -0.891  0.37310    
IslandDS       -5.162e-01  1.865e-01  2.916e+03  -2.767  0.00569 ** 
IslandFR        1.871e-01  2.876e-01  2.198e+03   0.650  0.51549    
fYear          -2.900e-02  5.796e-03  1.679e+03  -5.004 6.22e-07 ***
SexMale         2.589e+00  4.624e-02  3.176e+03  55.991  < 2e-16 ***
IslandAR:fYear -4.495e-02  1.871e-02  3.681e+03  -2.403  0.01632 *  
IslandCE:fYear  2.268e-03  1.354e-02  4.598e+03   0.168  0.86696    
IslandDS:fYear  4.375e-02  1.350e-02  4.156e+03   3.241  0.00120 ** 
IslandFR:fYear  2.893e-02  4.134e-02  2.172e+03   0.700  0.48407    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) IslnAR IslnCE IslnDS IslnFR fYear  SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR    -0.123                                                               
IslandCE    -0.198  0.099                                                        
IslandDS    -0.315  0.101  0.190                                                 
IslandFR    -0.225  0.072  0.131  0.197                                          
fYear       -0.551  0.177  0.302  0.464  0.344                                   
SexMale     -0.118 -0.015  0.005 -0.033 -0.026 -0.015                            
IslndAR:fYr  0.115 -0.873 -0.083 -0.092 -0.070 -0.200  0.013                     
IslndCE:fYr  0.174 -0.088 -0.910 -0.152 -0.112 -0.299 -0.007  0.096              
IslndDS:fYr  0.218 -0.061 -0.121 -0.851 -0.195 -0.359  0.012  0.068  0.107       
IslndFR:fYr  0.074 -0.023 -0.043 -0.099 -0.857 -0.135  0.024  0.028  0.040  0.217
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0760292 (tol = 0.002, component 1)

#predicting the model
IWLlmerdata <- ggpredict(IWLlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,WingLength=predicted,Island=group,Sex=facet)
IWLlmerdata2 <- merge(IWLlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11)) %>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear)%>% summarise(WingLength = mean(WingLength, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))

#plotting the model
IWLmod <- ggplot(IWLlmerdata2, aes(x = BirthYear, y = WingLength, fill = Island)) + 
  geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=WingLength,colour=Island), size = 2) +
  geom_errorbar(data=ADD_mean,aes(ymin=WingLength-WLsd,ymax=WingLength+WLsd, colour=Island),alpha=0.7) +
  geom_smooth(aes(colour = Island),method = "lm", se = FALSE) + 
  geom_ribbon(aes(ymin=WingLength-std.error,ymax=WingLength+std.error),alpha=0.5) +
  labs(x = "Birth Year", y= "Wing Length (mm)") +
  xlim(1990,2022) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind() + scale_fill_colorblind()

Head + Bill Length

#### Head + Bill ----
IHBlmer <- lmer(HeadBill ~  Island*fYear + Sex  + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IHBlmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: HeadBill ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +  
    (1 | BirdID)
   Data: ADD_I

REML criterion at convergence: 9601.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1841 -0.4234  0.0279  0.4642  5.4733 

Random effects:
 Groups   Name        Variance Std.Dev.
 BirdID   (Intercept) 0.22985  0.4794  
 Observer (Intercept) 0.04913  0.2217  
 Ageclass (Intercept) 0.04152  0.2038  
 Residual             0.14507  0.3809  
Number of obs: 5937, groups:  BirdID, 3120; Observer, 80; Ageclass, 5

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     3.588e+01  1.085e-01  7.815e+00 330.737   <2e-16 ***
IslandAR        4.153e-01  2.063e-01  1.204e+03   2.013   0.0443 *  
IslandCE       -1.041e-01  9.292e-02  3.741e+03  -1.121   0.2624    
IslandDS        8.141e-02  7.652e-02  3.252e+03   1.064   0.2874    
IslandFR       -7.000e-02  1.160e-01  2.381e+03  -0.603   0.5464    
fYear          -1.237e-03  2.492e-03  2.573e+03  -0.496   0.6196    
SexMale         8.420e-01  2.075e-02  3.016e+03  40.579   <2e-16 ***
IslandAR:fYear -2.977e-02  1.559e-02  5.975e+02  -1.910   0.0566 .  
IslandCE:fYear -8.098e-03  5.841e-03  3.978e+03  -1.386   0.1657    
IslandDS:fYear  1.404e-03  5.456e-03  3.918e+03   0.257   0.7969    
IslandFR:fYear  1.196e-02  1.638e-02  1.925e+03   0.730   0.4655    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) IslnAR IslnCE IslnDS IslnFR fYear  SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR    -0.066                                                               
IslandCE    -0.177  0.079                                                        
IslandDS    -0.277  0.088  0.233                                                 
IslandFR    -0.197  0.067  0.163  0.243                                          
fYear       -0.444  0.155  0.351  0.525  0.386                                   
SexMale     -0.089 -0.017 -0.001 -0.036 -0.029 -0.021                            
IslndAR:fYr  0.046 -0.944 -0.049 -0.062 -0.051 -0.133  0.009                     
IslndCE:fYr  0.152 -0.064 -0.918 -0.188 -0.138 -0.333 -0.002  0.047              
IslndDS:fYr  0.187 -0.049 -0.147 -0.851 -0.220 -0.393  0.012  0.041  0.131       
IslndFR:fYr  0.062 -0.021 -0.054 -0.113 -0.848 -0.148  0.023  0.019  0.050  0.218

#predicting the model
IHBlmerdata <- ggpredict(IHBlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,HeadBill=predicted,Island=group,SexEstimate=facet)
IHBlmerdata2 <- merge(IHBlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 22)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear)%>% group_by(Island,BirthYear)%>% summarise(HeadBill = mean(HeadBill, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))

#plotting the model
IHBmod <- ggplot(IHBlmerdata2, aes(x = BirthYear, y = HeadBill, fill = Island)) + 
  geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=HeadBill,colour=Island), size = 2) +
  geom_errorbar(data=ADD_mean,aes(ymin=HeadBill-HBsd,ymax=HeadBill+HBsd, colour=Island),alpha=0.7) +
  geom_smooth(aes(colour = Island),method = "lm", se = FALSE) + 
  geom_ribbon(aes(ymin=HeadBill-std.error,ymax=HeadBill+std.error),alpha=0.25) +
  labs(x = "Birth Year", y= "Head + Bill Length (mm)") +
  xlim(1990,2022) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind() + scale_fill_colorblind()

Arrange and plot

Imods <- ggarrange(IRTmod,IBMmod,IWLmod,IHBmod, common.legend=TRUE, labels = c("A","B","C","D"))
Imods

Lifespan and Reproduction

Preparing data

# make a non-redundant data set for lifespan and reproduction
ADD_Real <- ADD_I %>% filter(Island == "CN") %>% group_by(BirdID) %>% mutate(across(c(RightTarsus,BodyMass,WingLength,HeadBill),mean, na.rm=TRUE)) %>% filter(!duplicated(BirdID)) %>% select(RightTarsus,BodyMass,WingLength,HeadBill,BirdID,SexEstimate,BirthYear,Lifespan,SurvStatus,ReproductiveOutput) %>% mutate(SexEstimate = as.factor(SexEstimate))
# filter for missing values
ADD_tidy <- ADD_Real %>% filter(!is.na(RightTarsus) & !is.na(BodyMass)  & !is.na(WingLength)  & !is.na(HeadBill) & !is.na(SexEstimate)  & !is.na(Lifespan) ) %>% mutate(BodyCondition = BodyMass/RightTarsus) %>% mutate(SexEstimate = as.factor(SexEstimate))

Survival analysis with coxme

Model comparisons

###Lifespan analysis with coxme----

lscox0 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ BodyMass + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox1 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox2 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus*poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox3 <- coxme(Surv(Lifespan,SurvStatus) ~ poly(BodyCondition,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox4 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2)*SexEstimate + WingLength + HeadBill + (1|BirthYear), data= ADD_tidy)

AIC(lscox0,lscox1,lscox2,lscox3,lscox4) # lscox4 best

BIC(lscox0,lscox1,lscox2,lscox3,lscox4) # lscox4 best

# lscox4 was selected based on the best AIC and BIC values

lscox4

Cox mixed-effects model fit by maximum likelihood
  Data: ADD_tidy
  events, n = 1366, 1366
  Iterations= 10 75 
                   NULL Integrated    Fitted
Log-likelihood -8500.56  -8298.841 -8240.677

                   Chisq    df p    AIC    BIC
Integrated loglik 403.44  9.00 0 385.44 338.46
 Penalized loglik 519.77 35.45 0 448.86 263.81

Model:  Surv(Lifespan, SurvStatus) ~ RightTarsus + poly(BodyMass, 2) *      SexEstimate + WingLength + HeadBill + (1 | BirthYear) 
Fixed coefficients
                                        coef    exp(coef)   se(coef)     z       p
RightTarsus                       0.01881809 1.018996e+00 0.04997409  0.38 7.1e-01
poly(BodyMass, 2)1              -10.05540461 4.295298e-05 2.88625418 -3.48 4.9e-04
poly(BodyMass, 2)2                5.35344959 2.113361e+02 1.84019707  2.91 3.6e-03
SexEstimate1                      0.54043577 1.716755e+00 0.10887255  4.96 6.9e-07
WingLength                        0.02294382 1.023209e+00 0.02277508  1.01 3.1e-01
HeadBill                          0.00173566 1.001737e+00 0.05698445  0.03 9.8e-01
poly(BodyMass, 2)1:SexEstimate1 -13.89794735 9.208696e-07 3.49003613 -3.98 6.8e-05
poly(BodyMass, 2)2:SexEstimate1   2.33161533 1.029456e+01 2.56984891  0.91 3.6e-01

Random effects
 Group     Variable  Std Dev   Variance 
 BirthYear Intercept 0.7408611 0.5488751

Predicting and plotting coxme

####predict lifespan----
ADD_round <- ADD_Real %>%mutate(across(c(RightTarsus, BodyMass, HeadBill, WingLength), round, 0)) %>% distinct(RightTarsus, BodyMass, HeadBill, WingLength,SexEstimate) %>% filter(!is.na(RightTarsus) & !is.na(BodyMass)  & !is.na(WingLength)  & !is.na(HeadBill) & !is.na(SexEstimate))
ADD_end <- bind_rows(ADD_tidy,ADD_round)
rs <- predict_coxme(lscox4,newdata = ADD_end,type="risk", se.fit=TRUE)
rsend <- ADD_end %>% add_columns(rr = rs$fit[,1], se.fit = rs$se.fit[,1])
ADD_tail <- rsend %>% tail(nrow(ADD_round)) %>% mutate(SexEstimate = as.factor(SexEstimate)) %>% group_by(BodyMass,SexEstimate) %>% summarise(mean_rr=mean(rr),mean_se.fit=mean(se.fit))

# plotting coxme
coxmeplot <- ggplot(ADD_tail,aes(x=BodyMass,y=mean_rr)) + 
  geom_line(aes(colour=SexEstimate)) +
  geom_ribbon(aes(fill=SexEstimate,ymin=mean_rr-mean_se.fit,ymax=mean_rr+mean_se.fit), alpha = 0.25) + 
  labs(x = "Body Mass (g)", y = "Mortality Risk Score") +
  guides(fill=guide_legend(title="Sex"), colour=guide_legend(title="Sex")) +
  theme_tufte(base_size = 15, base_family = "Arial") + 
  scale_color_colorblind(labels = c("Female", "Male")) +
  scale_fill_colorblind(labels = c("Female", "Male"))

Reproductive Output using GLM

Model comparisons

###Reproduction glm----
ROglmF <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill  ,data = ADD_tidy, family = "poisson")
ROglmR <- glm(ReproductiveOutput ~ poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm0 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill  + Lifespan + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm1 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength  + HeadBill  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm2 <- glm(ReproductiveOutput ~ poly(RightTarsus,2)+ BodyMass + WingLength + HeadBill  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm3 <- glm(ReproductiveOutput ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm4 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + poly(WingLength,2) + HeadBill  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm5 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + poly(HeadBill,2)  + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")

AIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # ROglm4 slightly better than ROglm1 second

BIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # ROglm1 best 

#ROglm1 was selected based on AIC and BIC scores
vif(ROglm1)

                      GVIF Df GVIF^(1/(2*Df))
RightTarsus       3.252977  1        1.803601
BodyMass          2.553710  1        1.598033
WingLength        2.550319  1        1.596972
HeadBill          2.274822  1        1.508251
poly(Lifespan, 2) 1.212419  2        1.049333
SexEstimate       3.638709  1        1.907540
BirthYear         1.266770  1        1.125509

summary(ROglm1)

Call:
glm(formula = ReproductiveOutput ~ RightTarsus + BodyMass + WingLength + 
    HeadBill + poly(Lifespan, 2) + SexEstimate + BirthYear, family = "poisson", 
    data = ADD_tidy)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)        55.943199   6.971746   8.024 1.02e-15 ***
RightTarsus        -0.004860   0.032759  -0.148  0.88207    
BodyMass            0.075016   0.028836   2.602  0.00928 ** 
WingLength          0.015426   0.015372   1.004  0.31561    
HeadBill           -0.002388   0.038924  -0.061  0.95108    
poly(Lifespan, 2)1 23.108050   1.113788  20.747  < 2e-16 ***
poly(Lifespan, 2)2 -5.454150   0.702063  -7.769 7.93e-15 ***
SexEstimate1       -0.042798   0.066926  -0.639  0.52251    
BirthYear          -0.028461   0.003275  -8.692  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2018.43  on 736  degrees of freedom
Residual deviance:  879.58  on 728  degrees of freedom
  (629 observations deleted due to missingness)
AIC: 3146.9

Number of Fisher Scoring iterations: 5

Predicting and plotting GLM

####predict reproduction----
# grouping variables for glm plot
ROplotdata <- ADD_tidy %>% mutate(BodyMass = case_when(BodyMass > 16.5 ~ "16.5", BodyMass < 14.5 ~ "14.5",BodyMass < 16.5 & BodyMass > 14.5 ~ "15.5")) %>% filter(!is.na(BodyMass))

# predicting glm 
ROglmdat <- ggpredict(ROglm1, terms = c("Lifespan","BodyMass")) %>% rename(Lifespan=x,ReproductiveOutput=predicted,BodyMass=group)

#plotting glm
ROglmplot <- ggplot(ROglmdat, aes(x=Lifespan,y=ReproductiveOutput)) + 
  geom_point(data = ROplotdata, mapping=aes(x=Lifespan,y=ReproductiveOutput,colour=BodyMass),position="jitter") +
  geom_smooth(aes(colour=BodyMass), se = FALSE) + 
  geom_ribbon(aes(ymin=conf.low,ymax=conf.high,fill=BodyMass), alpha = 0.25) +
  labs(x = "Lifespan (years)", y= "Total Number of Offspring") +
  theme_tufte(base_size = 15, base_family = "Arial") +
  scale_color_colorblind() + scale_fill_colorblind()

Arrange and plot

fitplot <- ggarrange(coxmeplot,ROglmplot,  labels = c("A","B"))
fitplot

Time taken to run the whole script

Time difference of 38.09 secs

Acknowledgment

Thank you for going through the code and supporting open research. This paper would not have been possible without all authors contributing towards it. Massive thanks to everyone that has shared their excitement for this paper.